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Abstract. We investigate the interplay of diffraction and nonlinear effects during 
propagation of very short light pulses. Adapting the factorization approach to the 
problem at hand by keeping the transverse-derivative terms apart from the residual 
nonlinear contributions we derive an unidirectional propagation equation valid for 
weak dispersion and reducing to the slowly-evolving-wave-approximation for the case of 
paraxial rays. Comparison of numerical simulation results for the two equations shows 
pronounced differences when self-focusing plays important role. We devote special 
attention to modelling propagation of ultrashort terahertz pulses taking into account 
diffraction as well as Kerr type and second order nonlinearities. Comparing measured 
and simulated wave forms we deduce the value of the nonlinear refractive index of 
lithium niobate in the terahertz region to be three orders of magnitude larger than in 
the visible. 
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1. Introduction 

Improving theoretical description of propagation of light pulses consisting of just a few 
cycles has received considerable interest in the past two decades [D [21 El n El H m El 
lEniE]. However, relatively small number of investigations consider the variation in 
the transverse directions i.e. the diffraction effects. For example, analytical solutions 
are presented in [3] for diffraction and dispersion effects but only for the case of linear 
propagation. Our aim is to consider both nonlinear effects and transverse variation 
leading to broadening in the transverse direction but also self focusing for high intensities 
for the case when dispersion effects including attenuation have a minor role. 

Introducing additional dependence on one or both transverse coordinates makes 
mathematical modelling and numerical simulations more involved but is necessary for 
studying effects like self-guiding and hlament formation when competition between 
diffraction and nonlinearity plays important role |l2lll3lllHllHlll6] . Nonlinear envelope 
equations of the type derived in [1] were used for numerical simulation of shape effects 
[12] on the propagation of femtosecond pulses in media with weak dispersion and for 
studying hlament stability [H]. Analytical approach using soliton solutions can also 
be used for hlament propagation description [16] and collision-dynamics study of otical 
pulses [15] for certain hxed values of the nonlinear coefficient. Taking into account both 
the transverse (as well as longitudinal) prohle and nonlinearity is essential in describing 
interaction of optical and matter-wave soliton solutions characterizing trapped-atom 
arrangements PI- 

Examination of the term containing the transverse derivatives reveals that in 
the traditional slowly-varying-envelope approximation an implicit expansion in the 
transverse components of the wave vector is made by assuming that they are much 
smaller than the longitudinal component pointing in the propagation direction PI- 
This approximate treatment of transverse dependence also characterizes for example 
the slowly-evolving-wave approximation (SEWA) introduced in [1] as pointed out in 
[6]. We avoid this expansion pertaining to paraxial-ray approximation by modifying 
the factorization approach of Kinsler [10] to obtain an evolution equation suitable for 
numerical calculation which properly takes into account dihraction including space-time 
focusing ehects. 

2. Theoretical considerations 

To obtain an unidirectional propagation equation we follow the approach outlined in 
[To] and discussed earlier in more detail in [18]. We represent the physical electric 
held Eph{r,t) by analytic complex held E{r,t) whose Fourier transform contains only 
positive frequencies [191120] : 

Eph{r,t) = ^ [E{r,t) + c.c.]. 


( 1 ) 
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Taking the propagation direction to be the z axis the 3-diniensional wave eqation can 
be written as 

1 4-7r If 

{a'i + Vi) £:(r, i) - * E(r, t) = r, t) 

+ IMfIIffI + VV-E(r,t) ( 2 ) 

where -k denotes convolution, cl contains the linear part of the material response, 
Pnl {E,r,t) is the nonlinear polarization, J{r,t) the current density and we assumed 
a non-magnetic material. Transforming into wave-vector and frequency space we can 
write ([2]) as 


{-kl-kl + f) E{k,i^) = -Q, 


(3) 


with —Q denoting all terms corresponding to the right-side of ([2]) and ff^ = uj‘^eL{oj)/c^. 
Note that in contrast to [10] we do not place the term with transverse derivatives in the 
residual term Q but keep it exlicitly. Performing the factorization in ([3]): 


K - - k]^ (fc, + - k]^ Eik,u) = Q, 


we can express as a sum of two terms: 
E{k,u!) = 


- k]_ \K - - k]_ kz + 13'^ - k] 


Dehning E^^'> by 


E^^\k,cj) = 


±1 


-Q 


we can write E = E^^'> -1- E^~'^ where E^^^ obey the equations: 


k,TJ^-kl)E^^\k,uj) = 


±1 


-Q. 


Q 


(4) 


(5) 


( 6 ) 


(7) 


Transforming back to ordinary space with respect to kz (but not k±) we obtain the two 
equations governing propagation of k±,u): 


dzE^^\z,k±,u) = Ti-vZ/d^ — kfE^'^\z,k±,uj) ± 


2v^ 


kl 


Q 


( 8 ) 


Equations (jS]) can be regarded as an alternative derivation of the unidirectional pulse 
propagation equation of Kolesik and Moloney which is based on modal expansion of 
the helds. A potential disadvantage of relying on result ([H]) is that one has to work in 
the perpendicular wave vector space and without expansion one can not meaningfully 
transform to perpendicular derivatives in the ordinary space. Since our hnal intention is 
to arrive at an unidirectional propagation equation using only the forward propagating 
held we have to mention an important constraint arising from attenuation taken 
into account through the (positive) imaginary part of 13“^ {oj). That term provides 
suppression for the forward propagating component but as seen form ([8]) it enhances the 
backward propagating term. Poor phase matching is expected to provide suppression 
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for backward propagation but for long enough propagation distances compared to 
attenuation length that component may be amplihed enough in order to make its 
effect through nonlinearity non-negligible. As a consequence, for the applicability of 
the unidirectional equation we require the condition of weak dispersion to be satished 
for the full frequency spectrum of the pulse and the propagation distance to be shorter 
than the attenuation length. This is in accordance with constraints on the applicability 
of the slowly-evolving-wave-approximation |T] concerning weak dispersion and restricted 
propagation length as discussed in [21] and section 11.6 of [22] . 

For a general transverse dependence one can use the cartesian components of the 
transverse wave vector kx, ky and the corresponding eigenfunctions eyrp{ikxX+ikyy) as in 
[7] . However, since our main interest resides in cylindrically symmetric pulse prohles we 
exploit that symmetry to make the numerical modelling less computationally intensive 
and use the simulation with arbitrary transverse dependence only as a check of numerical 
procedures. 

We make the transformation to transverse wave vector space using the fact that 
the solution of the eigenvalue problem 

Vl/(r,0) + fcl/(r,0) = O, (9) 


with r, 0, z being the cylindrical coordinates, is 

OO 

/(r, + B.„ iV„(nr)] e-’"^ 


( 10 ) 


m=0 


where Jm is the Bessel function and Nm is the Neumann function of order m. Assuming 
uniform medium in the transverse direction and regularity of the solution at r = 0 with 
rotational symmetry around the z axis means that only the Jo term survives. We can 
thus expand the transverse dependence as 


/(r) = / A(fcx) Mk±r)dk^ 

Jo 

and using the integral relation 

/ rJo{kfr)Jo{k^r)dr = —6{kf-kx), 

Jo 

where S{k) is the Dirac delta function, we can invert (ITTil to get 

poo 

A{k±) = k± rf{r)Jo{k±r)dr. 


( 11 ) 


( 12 ) 


(13) 


There is an alternative approach to expression oa which relies on explicit 
compactihcation of the transverse coordinate. Even without a strict boundary condition 
that the electric held vanishes at a hxed r value, we may require that it be zero at some 
(large) value of r = a. That condition restricts the continuous integration variable in 
(na to discreet values kf^ satisfying the condition Jo^kfa) = 0 and thus leads to 

CXD 

fir) = '^AmMkir). 

m=l 


( 14 ) 
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Denoting the zeros of the Bessel function Jo{x) by rn 
kfn = Xm/cL- Usiug the orthogonality relation [2^ 

^2 

rJo{xmr/a) Jo{xnr/a)dr = 
one can compute the expansion coefficients: 



1,2,3,... we can write 


(15) 


2 r 

Am = -iT-rt —t? / rf{r)Jo{xmr/a)dr. (16) 

a^Ji[XmY Jo 

We used this approach as a check of numerical implementation based on oa and 
indeed ascertained that in our numerical simulations the above two approaches give 
indistinguishable results if the value a of the radius at which the field vanishes is chosen 
large enough and for considered limited propagation distances. We acknowledge that the 
above compactification should be regarded only as an approximate numerical procedure 
and used with care since it leads to violation of general integral properties [25l |26l [2^ 
of angular-spectrum representation for diffracted wave fields. 

If the residual term on the right side of ([7]) does not mix the forward propagating 
field EA) with the backward propagating E^A considerably, i.e. it is small compared 
to propagation terms not mixing the two field components then one can decouple the 
corresponding equations for forward and backward propagation. 

When taking into account both forward and backward propagating fields 
introduction of envelopes with slower spatial (in the 2 ; direction) variation does not 
bring any real advantage since that introduces fast varying contributions in the residual 
term which mixes the two fields. 

If the mixing of forward and backward propagating fields through the residual term 
is small enough so that a unidirectional propagation applies to good accuracy it is 
advantageous to introduce the envelope A{z,r_i,t) with slower 2 ; and t variation than 
the (complex) field itself through 


^ A{z,rs_,t) exp[i(/?o2 - wo^)], 


(17) 


where ojo is the predetermined central (carrier) frequency and (Jq = n{uJo)oJo/c. In this 
way dH]) for the forward propagating field is replaced by the following expression: 

dzA{z,k^,u) = -k\- IJo)M^^ 

+ — , ^ Qa{A), (18) 

2^/32 - kl ^ ^ ^ ^ 

with the residual term on the right side depending also on the envelope. It is customary 
to introduce the moving frame by transforming to new time r and longitudinal distance 
C variables [1]: 


r = t- Yiz, C = z, 


(19) 


with Vg{uo) = l//3i the group velocity at Uq. The unidirectional envelope evolution 
equation (USD then becomes 


d(;A{C,k±,u) = i 




kl 


(Jo — (Jl^ 


A{C,k±,u) 
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Qa(A). 


( 20 ) 


3. Numerical simulation 


We now turn to numerical simulation of light pulse propagation. As a first example we 
consider the case examined in [1] but without diffraction effects being taken into account 
there. A Kerr type nonlinearity is assumed with electric held and nonlinear polarization 
having the same direction. This is a reasonable approximation since the nonlinear effect 
in itself is small and including a small correction in the form of projection on the electric 
held vector is not expected to be signihcant In our simulations we do distinguish 
the transverse and longitudinal components of the electric held with respect to the 
axis of propagation based on values of the length of wave vector and its transverse 
component and add up these two components of electric held separately. However, 
for the cases studied numerically in the following we observe that this separation of 
transverse and longitudinal components gives noticeable diherence compared to the case 
when neglecting it only for large enough transverse-coordinate values where typically 
the electric held magnitude drops to 1-2% of its central value. After transforming to u 
(i.e. fdioj)) and /c^ space the scalar form of fl20|) is thus appropriate: 


d(^A{C,k^,u:) = i 


-kf - /do- fdiuj 


A(C,fc^,n;) 


2m{uj ujoY 


B/A), 


c‘^\//d‘^ — kf 

where the nonlinear polarization has been written as 


( 21 ) 


Pnl{z, r_L, t) = Anl{z, r_L, t, A) exp[i(/3o^ - Wot)] + c.c., 

and B{A) is the Fourier transform of Ajsfi to frequency and perpendicular wave vector 
space. If we expand f/f^ — kf in the hrst term of fl2Tl) to linear order in kf and in the 
term containing kf approximate the refractive index n{u) with its value at the carrier 
frequency uo we recover the linear terms of the SEWA equation (6) in [1]. We also recover 
the nonlinear term in that equation if in our nonlinear term we put k± to zero and again 
take the refractive index frequency independent and equal to its value at uq. Simulation 
results show that these simplihcations are quite good approximations for the case 
considered in jl] consisting of a pulse with central wavelength Aq = 0.8 /im propagating 
in fused silica with a hyperbolic-secant shaped envelope As(t) = 1/ cosh [1.76f/rp] and 
Tp = 2.67 fs. However, pronounced difference characterizes the solutions of the two 
equations when nonlinearity induced self focusing plays important role. In order to 
allow for dispersion and diffraction effects to develop over larger propagation distance 
we decreased the peak intensity used in [1] by a factor of two to 2 x 10^^ W/cm^ 
for analyzing propagation in a Kerr medium with Ajvl = (27r)“^non2|ApA where 
n2 = 3 X 10“^^ cm^/W is the nonlinear index of refraction and |Ap is normalized 
to give the intensity. 
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We observe that considering the propagation to be lossless for the above pulse 
properties and propagation distances not exceeding 100 /rm is a very good approximation 
since the imaginary part of the refractive index of considered silica glass in the frequency 
interval where the spectral distribution of the pulse is not less than 2% of its maximum 
is not exceeding 0.0004 |2H] implying an attenuation length larger than 1500 pm and 
even at the cut-off wavelength Xm = 5.3 pm (see caption of hgure [T]) it is 140 pm. 
In case of strong dispersion or propagation distances exceeding the attenuation length 
more careful treatment is required prefering solution of the full nonlinear wave equation 

m- 

We start by analyzing propagation of a beam with Gaussian transverse distribution 
whose amplitude is proportional to exp[—r^/(2a^)] with = 3 pm. The transverse- 
derivative term in expression (6) of [T] contains the term uj + ujo in the denominator and 
thus diverges when the envelope spectrum extends close to —cjq which can be the case for 
single-cycle pulses. It is also present in the studied example which means that in order 
to avoid overflow during iterations one has to implement a cut-off for u + uq values 
approaching zero. Our simulations show that this value can be signihcantly smaller 
than ujo (for example O.Icjq) without causing numerical problems. In order to avoid 
underestimation of the electric held the cut-off should be signihcantly smaller than the 
value at which the envelope falls to one half of its maximum as illustrated by hgure [H In 
that hgure we compare results for diherent cut-oh values after propagation distances of 
40 pm and 80 pm at central (r = 0) position. We observe similar relative diherences for 
oh-axis positions. For a stronger focused beam with Or = 1 pm the diherences are even 
larger. In the following comparisons we use the cut-oh value = 5.3 pm which leaves 
out only 0.5% of the envelope spectral distribution. We remark that in our result 
the propagation condition fl > k± applies and in the diverging nonlinear term there is 
no need to impose constraint on the lower limit of \/3 — k±\ due to the integrable nature 
of the singularity. 

Next, we compare the transverse prohles of the pulse by plotting the huence 
as function of transverse coordinate after diherent propagation distances using the 
solution of fl2ll) and of the corresponding expression (6) of [T]. We calculate the 
huence by integrating \E{z,r,u})\‘^ over the frequency u. Figure [2] shows the transverse 
huence prohles for diherent propagation distances starting with a Gaussian amplitude 
distribution exp[—r^/(2a^)] with = 3.5 pm. In hgure[2]we see the ehect of self focusing 
which is dominant if the intensity is sufficiently high and the focusing of incoming pulse 
not too strong. The inhuence of the initial shape of the pulse on the propagation is 
quite important similarly to the case when femtosecond pulses propagate in air nasi]. 
We remark that in case of studied pronounced self focusing the resultant huence can 
be quite high and not far from but not reaching the damage limit for fused silica |30j . 
In hgure [3] we compare the time dependence of the electric held of the pulse after 
propagating 62 pm at r = 0 using the two approaches. Our results are in general close 
to the results obtained by using the SEWA approach. However, in case of pronounced 
self focusing corresponding to propagating z ^ 60 pm we observe that SEWA leads 
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Figure 1. Electric field at z = 40 ym (upper panel) and z = 80 ym (lower panel) 
for r = 0. The full line corresponds to the cut-off wavelength Am = 5.3 ym while the 
dash line is for Am = 1-04 ym which corresponds to the (lower) frequency at which the 
envelope takes one half of its maximum value. 



Figure 2. Transverse fluence profiles for initial distribution with Ur = 3.5 ym. Our 
result (a) and using the result from integrating the corresponding equation (6) of [T] 
(b). 


up to 20% larger central fluence values compared to our result. It is not difficult to 
understand the cause of this behaviour by recalling that SEWA amounts to dropping 
the transverse component of the wave vector in the square root of the nonlinear term in 
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Figure 3. Time dependence of the electric field at r = 0 after propagation of 62 gm. 
Solid line is our result and dash line is obtained by using the result of the equation 
from [T]. 


equation fl20|) thus diminishing the magnitude of that contribution for held components 
with nonzero transverse coordinate. This in turn leads to relative enhancement of the 
nonlinear effect for the central region of pulse. We can see from hgure |3] that this 
enhancement comes from larger electric-held amplitude in the descending part of the 
envelope while the self-phase-modulation ehect does not show noticeable diherence. We 
also observe pronounced diherence between the two approaches for large enough values 
of the transverse coordinate where the electric held is typically two orders of magnitude 
smaller than at central position. In this oh-axis region SEWA typically underestimates 
the electric-held magnitude. 

We now turn to analyzing terahertz single-cycle pulse propagation as exemplihed by 
observation of nonlinear lattice response in |3T] where a pulse with peak intensity of 100 
MW/cm^ was transmitted through a 2 mm thick LiNbOs crystal cooled to 80K. Based 
on the measured incoming (reference) electric held strength we construct a complex 
envelope function by retaining only the positive-frequency part of the Fourier transform 
and then introducing a shift by the carrier frequency Uq = 3.75 ps“^ obtained using 
the formula Uq = u\E{u)\‘^du/ \E{u)\‘^du. The measured reference electric held 

and the one calculated using the envelope function, together with its shape are shown 
in hgure 01 In hgure [5] the observed transmitted held strength is shown next to the 
reference one and the simulated result assuming linear propagation with and without 
absorption and with frequency-dependent refractive index of the form n = A+b f‘^+C f^ 
with A = 4.73, B = 1.510“^, C = 8.510“^° and / = u:/{2ttc) with ca/c expressed in 
cm“^ [32] . 

In this example absorption is not completely negligible but has a small ehect as 
seen from hgure [5] and measurements reported in [32]. Comparing the incoming and the 
transmitted pulse we conclude that the absorption coefficient at the carrier frequency 
should be a{uJo) ~ 1 cm“^ and taking into account the measured increase with frequency 
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Figure 4. The envelope function (dash line) and the measured electric field (solid 
line) together with the electric field calculated using the envelope function. 



Figure 5. The measured original (dash-dot line), the measured transmitted electric 
field (solid line) and the calculated transmitted field assuming linear propagation with 
dispersion from [32] without absorption (dot line) and with absorption coefficient 
a = 1 cm“^ (dash line). 


[32] we arrive at the value a{um) <1.5 cm“^ where Um ~ 8 ps“^ is the frequency at 
the effective upper limit of the frequency spectrum of the pulse. Since the real part 
of the propagation constant fir{^m) ~ 1200 cm“^ we ascertain that the condition of 
weak dispersion is satisfied at the high end of spectrum and also for smaller frequencies 
characterizing the pulse (measurements extending down to 0.25 THz did not show 
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increase of absorption coefficient [53] )• We also observe that the attenuation length 
is at least three times the propagation distance for pulse frequencies. 

Introducing a third order Kerr type nonlinearity clearly improves the agreement 
of the simulated and measured results, notably the increased distance between the 
hrst positive and hrst negative lobes as can be seen in hgure El In this case we did 
not consider diffraction effects by assuming a uniform transverse distribution. The 
value of the nonlinear index of refraction giving the closest overall agreement with 
measurement corresponds to ^2 ~ 4 x 10“^^ cm^/W which is almost four orders of 
magnitude larger than its value in the visible part of spectrum. We also included the 
second-order nonlinear effect giving rise to the complex nonlinear-polarization envelope 



A^^{{z,r±,t) = ^{A{z,r±,tY exp[i{f3oz - ujot)] 
Zn 

+ 2|/1(2, rx.OP exp|-i(/3o2 - n'o*)]} 


( 22 ) 


with the restriction that only the positive-frequency part of the polarization should be 


taken by enforcing the constraint u > —ojq for the envelope. For the nonlinear optical 
coefficient ^33 we use its relationship to the near-infrared refractive index and clamped 
electro-optic coefficient r 33 (expression ( 6 ) of [35]) leading to 0^33 ~ —160 pm/V with 
^"33 = 30.8 [36]. The effect of this term is quite small for the considered intensity and 
propagation distance but nevertheless we observe that it brings the simulated electric- 
held variation slightly closer to the measured signal in the descending region of intensity. 



2 


3 

Time [ps] 


4 


5 


Figure 6. Comparison of transmitted and simulated results for different values of 
nonlinear index of refraction assuming uniform transverse distribution. 
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In order to model diffraction effects we take the incoming pulse with a Gaussian 
transverse prohle and for the measured values of reference and transmitted held take 
an average over the width of diameter D ^ 1 mm corresponding to the experimental 
condition. Using the estimated value of the width parameter ~ 0.6 mm shows 
marked improvement of agreement with measurement for the hrst oscillation and brings 
the preferred value of the nonlinear index of refraction to around 10“^^ cm^/W as shown 
in hgure [71 One has to acknowledge that the relatively small size of the nonlinear effect 
for the studied intensity and propagation distance introduces considerable uncertainty 
in that estimate and for more precise determination dedicated experiments with 
larger intensity and/or propagation distance would be required. In case of longer 
propagation distances surpassing the attenuation length a more careful consideration 
of attenuative dispersion would be required as shown by results obtained in [211 1^ - 
We also investigated the inhuence of frequency-dependent absorption coefficient slowly 



Time [ps] 


Figure 7. Comparison of transmitted and simulated results for different values of 
nonlinear index of refraction for Gaussian transverse profile of incoming beam with 
Ur = 0.6 mm. 

increasing with frequency motivated by results of [32] but since it did not produce 
noticeable change in the agreement with observation we only show results for the 
constant value of a = 1 cm“^. 


4. Conclusions 

We studied diffraction effects during nonlinear propagation of few-cycle light pulses 
with axial symmetry in the framework of unidirectional propagation equation derived 
for the case of weak dispersion by suitably modihed factorization method not relying on 
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the paraxial approximation. The slowly-evolving-wave approximation is obtained as a 
special case and numerical simulations show similar results with signihcant differences 
present in case of high intensity accompanied with pronounced self focusing effects. 
Analysis of short terahertz pulse propagation in LiNbOa shows both nonlinear and 
diffraction effects despite of short propagation distance and indicates presence of a Kerr- 
type nonlinearity with nonlinear index of refraction n 2 ~ 10“^^ cm^/W which is three 
orders of magnitude larger than its value in the visible part of spectrum. This points 
to signihcant contribution of lattice vibration anharmonicity and indicates pronounced 
suitability of terahertz pulses for investigating nonlinear lattice dynamics. 
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